home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / t3_1 / sources.lha / sources / sys / big_arith.t < prev    next >
Text File  |  1988-02-05  |  17KB  |  438 lines
  1. (herald big_arith
  2.   (env tsys (osys fixnum) bignum))
  3.  
  4. ;;; Addition:
  5.  
  6. ;;; Places the modular sum of u, v, and carry, in element i of
  7. ;;; vector destv.  Returns the new carry.
  8.  
  9. (define-integrable (add-step u v destv i carry)
  10.   (receive (sum carry)
  11.            (%digit-add u v carry)
  12.     (set (bignum-digit destv i) sum)
  13.     carry))
  14.  
  15. ;;; add-magnitudes takes two bignums, adds their magnitudes, and returns
  16. ;;; the magnitude of their sum.
  17.  
  18. (define (add-magnitudes u v)
  19.   (cond ((fx< (bignum-length u) (bignum-length v))
  20.          (really-add-magnitudes v u))
  21.         (else
  22.          (really-add-magnitudes u v))))
  23.  
  24. ;;; Assumes that first number is longer than second number.
  25.  
  26. (define (really-add-magnitudes u v)
  27.   (let* ((u-size (bignum-length u))
  28.          (v-size (bignum-length v))        ; u-size >= v-size
  29.          ;; The following test tries to predict whether there will be carryout.
  30.          (lc (if (or (%digit-greater? (bignum-digit u (fx- u-size 1))
  31.                                       *half-max-hyperdigit*)
  32.                      (and (fx= u-size v-size)
  33.                           (%digit-greater? (bignum-digit v (fx- v-size 1))
  34.                                            *half-max-hyperdigit*)))
  35.                  (fx+ 1 u-size)
  36.                  u-size))
  37.          (sum (create-bignum lc)))
  38.     (do ((carry 0 (add-step (bignum-digit u i) (bignum-digit v i) sum i carry))
  39.          (i 0 (fx+ 1 i)))
  40.         ((fx= i v-size)
  41.          ;; Propagate carry
  42.          (do ((carry carry
  43.                      (add-step (bignum-digit u i) 0 sum i carry))
  44.               (i i (fx+ 1 i)))
  45.              ((or (fx= i u-size) (fx= carry 0))
  46.               ;; Copy remaining digits from u
  47.               (do ((i i (fx+ 1 i)))
  48.                   ((fx= i u-size)
  49.                    (cond ((fx= 0 carry)
  50.                           (set-bignum-length! sum u-size))
  51.                          (else
  52.                           (set (bignum-digit sum i) 1)))
  53.                    sum)
  54.                 (set (bignum-digit sum i) (bignum-digit u i)))))))))
  55.  
  56. ;;; Subtraction:
  57.  
  58. (define-integrable (subtract-step u v destv i carry)
  59.   (receive (diff carry)
  60.            (%digit-subtract u v carry)
  61.     (set (bignum-digit destv i) diff)
  62.     carry))
  63.  
  64. ;;; Returns the difference of the magnitudes of two bignums.
  65. ;;; Assumes that first number is at least as long as than second.
  66.  
  67. (define (subtract-magnitudes u v)
  68.   (let* ((la (bignum-length u))
  69.          (lb (bignum-length v))
  70.          (diff (create-bignum la)))
  71.     (do ((carry 0 (subtract-step (bignum-digit u i)
  72.                                  (bignum-digit v i)
  73.                                  diff i carry))
  74.          (i 0 (fx+ 1 i)))
  75.         ((fx= i lb)
  76.          (do ((carry carry (subtract-step (bignum-digit u i) 0 diff i carry))
  77.               (i i (fx+ 1 i)))
  78.              ((or (fx= la i) (fx= carry 0))
  79.               (do ((i i (fx+ 1 i)))
  80.                   ((fx= la i)
  81.                    (bignum-trim! diff))
  82.                 (set (bignum-digit diff i) (bignum-digit u i)))))))))
  83.  
  84.  
  85. ;;; Multiplication:
  86.  
  87. ;;; Multiplies a bignum by a scalar, and adds the product into an
  88. ;;; accumulator beginning at the specified start index.
  89. ;;; Assumes that accum is at least of length v-size + start + 1.
  90.  
  91. (define (multiply-step scalar v v-size accum start)
  92.   (do ((carry 0 (receive (d1 d0)
  93.                          (%digit-multiply scalar (bignum-digit v j))
  94.                   (add-step (bignum-digit accum (fx+ 1 i))
  95.                             (fx+ d1 carry)  ; Note: d1 < *max-hyperdigit*
  96.                             accum
  97.                             (fx+ 1 i)
  98.                             (add-step (bignum-digit accum i) d0 accum i 0))))
  99.        (i start (fx+ 1 i))
  100.        (j 0 (fx+ 1 j)))
  101.       ((fx>= j v-size)
  102.        (cond ((fxn= carry 0)
  103.               (add-step (bignum-digit accum (fx+ 1 j))
  104.                         0 accum (fx+ 1 j) carry)))
  105.        accum)))
  106.  
  107. ;;; Similar in action to ADD-MAGNITUDES only with product.
  108.  
  109. (define (multiply-magnitudes u v)
  110.   (let* ((u-size (bignum-length u))
  111.          (v-size (bignum-length v))
  112.          (l-product (fx+ (fx+ u-size v-size) 1))       ;Why + 1?
  113.          (product (create-bignum l-product)))
  114.     (do ((i 0 (fx+ 1 i)))
  115.         ((fx= i u-size) (bignum-trim! product))
  116.       (multiply-step (bignum-digit u i) v v-size product i))))
  117.  
  118. ;;; Division:
  119.  
  120. ;;; The New Division Routine by J.R.
  121. ;;; (In a division u/v, u is the dividend and v is the divisor.)
  122.  
  123. (define (div2-magnitudes u v)
  124.   (let* ((m+n (bignum-length u))
  125.          (n (bignum-length v))
  126.          (d (fx- *bits-per-hyperdigit*
  127.                  (fixnum-howlong (bignum-digit v (fx- n 1)))))
  128.          (new-u (make-and-replace-bignum (fx+ m+n 1) u 0 0 m+n)))
  129.     ;; Normalize
  130.     (cond ((fx> d 0)
  131.            (really-div2-magnitudes (bignum-ashl! new-u d)
  132.                                    (bignum-ashl v d)
  133.                                    m+n n d))
  134.           (else
  135.            (really-div2-magnitudes new-u v m+n n 0)))))
  136.  
  137. ;;; Knuth equivalences:
  138. ;;;   u[j]     ->   (bignum-digit u m+n-j)
  139. ;;;   u[j+1]   ->   (bignum-digit u (- m+n-j 1))
  140. ;;;   u[j+n]   ->   (bignum-digit u (- m-j))
  141.  
  142. (define (really-div2-magnitudes u v m+n n d)
  143.   (let* ((v1 (bignum-digit v (fx- n 1)))
  144.          (v2 (if (fx> n 1) (bignum-digit v (fx- n 2)) 0))
  145.          (m (fx- m+n n))
  146.          (q (create-bignum (fx+ m 1))))
  147.     (do ((m-j   m   (fx- m-j   1))
  148.          (m+n-j m+n (fx- m+n-j 1)))
  149.         ((fx< m-j 0)
  150.          (return (bignum-trim! q)
  151.                  (bignum-trim! (bignum-ashr! u d))))
  152.  
  153.       ;; Calculate one digit of the quotient.
  154.       (let ((q^ (compute-q^ (bignum-digit u m+n-j)
  155.                             (if (fx> m+n-j 0) (bignum-digit u (fx- m+n-j 1)) 0)
  156.                             (if (fx> m+n-j 1) (bignum-digit u (fx- m+n-j 2)) 0)
  157.                             v1
  158.                             v2)))
  159.         (set (bignum-digit q m-j) q^)
  160.  
  161.         ;; Multiply and subtract: (uj...uj+n)  -:=  q^ * (v1...vn)
  162.         (iterate loop ((prev 0)
  163.                        (borrow 0)
  164.                        (i m-j)
  165.                        (k 0))
  166.           (cond ((fx< k n)
  167.                  (receive (d1 d0)
  168.                           (%digit-multiply q^ (bignum-digit v k))
  169.                    (receive (sum carry)
  170.                             (%digit-add d0 prev 0)
  171.                      (let ((d1 (fx+ d1 carry)))
  172.                        (receive (diff borrow)
  173.                                 (%digit-subtract (bignum-digit u i) sum borrow)
  174.                          (set (bignum-digit u i) diff)
  175.                          (loop d1 borrow (fx+ i 1) (fx+ k 1)))))))
  176.                 (else
  177.                  (receive (diff borrow)
  178.                           (%digit-subtract (bignum-digit u i) prev borrow)
  179.                    (set (bignum-digit u i) diff)
  180.                    (cond ((fxn= borrow 0)
  181.                           (add-back u v q q^ n m-j))
  182.                          ((fxn= diff 0)
  183.                           (error "inconsistency in bignum division!~%  ~S"
  184.                                  `(div ,u ,v)))
  185.                          )))))))))
  186.  
  187. ;;; In the loop, q^ starts out being 0, 1, or 2 larger than the actual
  188. ;;; first digit of the quotient of the bignums u and v.
  189. ;;; The loop may decrement q^, eliminating all cases where q^ is
  190. ;;; two larger, and most cases where it is one larger.
  191.  
  192. ;;; The initial guess for q^ is obtained by dividing u's first two digits
  193. ;;; by v's first digit.
  194. ;;; r^ is initially the remainder of the division, but as q^ is
  195. ;;; decremented, r^ maintains as its value the actual difference between
  196. ;;; the dividend  {u[j] u[j+1]}  and the product  q^*v[1].
  197. ;;; We stop pruning q^ as soon as its product with the second digit
  198. ;;; of v, exceeds r^ adjoined with the third digit of u.
  199.  
  200. (define (compute-q^ uj uj+1 uj+2 v1 v2)
  201.   (labels (((loop q^ r^)
  202.             ;; Test to see whether, in Knuth's notation,
  203.             ;;  q^*v[2] > r^*b + u[j+2]
  204.             ;;     where r^ = u[j]*b + u[j+1] - q^*v[1]
  205.             ;; We use the same tricks he does in his MIX code.
  206.             (receive (a1 a0)
  207.                      (%digit-multiply q^ v2)
  208.               ;; {a1 a0}  =  v[2]*q^
  209.               (cond ((and (not (%digit-greater? a1 r^))
  210.                           (or (fxn= a1 r^)
  211.                               (not (%digit-greater? a0 uj+2))))
  212.                      q^)
  213.                     (else
  214.                      (test (fx- q^ 1) r^)))))
  215.            ((test q^ r^)
  216.             ;; Adjust r^.  q^ must get no smaller if r^ overflows here.
  217.             (receive (sum carry)
  218.                      (%digit-add r^ v1 0)
  219.               (if (fxn= carry 0)
  220.                   q^
  221.                   (loop q^ sum)))))
  222.     ;; uj <= v1
  223.     (cond ((fx= uj v1)
  224.            ;; Note that in this case, uj+1 <= v2 also.  E.g. in
  225.            ;; (/ 8123,4567,... 8123,xxxx,...) -> q= FFFF, r= C68A (=8123+4567)
  226.            ;; we know that xxxx >= 4567.
  227.            (test *max-hyperdigit* uj+1))
  228.           (else
  229.            (receive (q^ r^)
  230.                     (%digit-divide uj uj+1 v1)
  231.              (loop q^ r^))))))
  232.  
  233. ;;; We come here if compute-q^ screwed up and gave us a bogus guess for
  234. ;;; the quotient digit.  The probability of this happening is about
  235. ;;;  2 / b  where b = (1+ *max-hyperdigit*).
  236.  
  237. (define (add-back u v q q^ n m-j)
  238.   (*let-us-know* u v q^)
  239.   (set (bignum-digit q m-j) (fx- q^ 1))
  240.   (iterate loop ((carry 0)
  241.                  (i m-j)
  242.                  (k 0))
  243.     (cond ((fx< k n)
  244.            (receive (sum carry)
  245.                     (%digit-add (bignum-digit u i) (bignum-digit v k) carry)
  246.              (set (bignum-digit u i) sum)
  247.              (loop carry (fx+ i 1) (fx+ k 1))))
  248.           (else
  249.            (set (bignum-digit u i) 0)))))
  250.  
  251. (lset *let-us-know*
  252.   (lambda (u v q^)
  253.     (if (experimental?)
  254.         (format (terminal-output)
  255.                 '("~&;Please send the following numbers to the T implementors.~%"
  256.                   ";This is not an error.~%; q^ = ~s~%; u  = ~s~%; v  = ~s~%")
  257.                 q^ u v))
  258.     (set *let-us-know* false)))
  259.  
  260. ;;; This depends on (>= *bits-per-hyperdigit* *bits-per-fixnum*)
  261. ;;; => T and, in %B-F-DIV2, on (FIXNUM-ABS most-negative-fixnum)
  262. ;;; => most-negative-fixnum which is true in T2.9 and T3.0 and should
  263. ;;; be true in general.
  264.  
  265. (define (b-f-multiply u v)
  266.   (let* ((u-size (bignum-length u))
  267.          (product (create-bignum (fx+ 1 u-size))))
  268.     (multiply-step (fixnum-abs v) u u-size product 0)
  269.     (set-bignum-sign! product
  270.                       (if (fx= (bignum-sign u)
  271.                                (if (fx< v 0) -1 1))
  272.                           1
  273.                           -1))
  274.     (normalize-integer (bignum-trim! product))))
  275.  
  276. ;;; Division of a bignum by a fixnum.
  277.  
  278. (define (b-f-div2 u v)
  279.   (cond ((fx= v 0)
  280.          (error "division by zero~%  (DIV ~s ~s)" u v))
  281.         (else
  282.          (%b-f-div2 u v t))))
  283.  
  284. ;;; An unnormalized result is required by the printer and INTEGER->POINTER.
  285.  
  286. (define-integrable (b-f-div2-unnormalized u v)
  287.   (%b-f-div2 u v nil))
  288.  
  289. (define (%b-f-div2 u v normalize?)
  290.   (let* ((v-abs (fixnum-abs v))
  291.          (d (fx- *bits-per-hyperdigit*
  292.                  (fixnum-howlong v-abs)))
  293.          (sign (if (fx= (bignum-sign u)
  294.                         (if (fx< v 0) -1 1))
  295.                    1
  296.                    -1)))
  297.     (cond ((fx= d 0)
  298.            (b-f-div2-aux u v-abs normalize? d sign))
  299.           (else
  300.            (b-f-div2-aux (bignum-ashl u d)
  301.                          (fixnum-ashl v-abs d)
  302.                          normalize?
  303.                          d
  304.                          sign)))))
  305.  
  306. (define (b-f-div2-aux u v normalize? d sign)
  307. ;  (format t "~&0> ~D ~D ~A ~D ~D~%" u v normalize? d sign)
  308.   (let* ((u-size (bignum-length u))
  309.          (q (create-bignum u-size)))
  310. ;    (format t "~&2> ~D ~D ~%" u-size q)
  311.     (iterate loop ((i (fx- u-size 1))
  312.                    (r 0))
  313.       (cond ((fx< i 0)
  314.              ;; Last time through, r is the remainder.
  315. ;             (format t "~&2> ~D ~%" q)
  316.              (set-bignum-sign! q sign)
  317. ;             (format t "~&3> ~D ~%" q)
  318.              (bignum-trim! q)
  319. ;             (format t "~&4> ~D ~%" q)
  320.              (return (if normalize? (normalize-integer q) q)
  321.                      (let ((r (fixnum-lshr r d)))
  322.                        (if (bignum-positive? u) r (- 0 r)))))
  323.             (else
  324.              (let ((word0 (bignum-digit u i)))
  325.                (receive (qi r)
  326.                         (%digit-divide r word0 v)
  327.                  (set (bignum-digit q i) qi)
  328.                  (loop (fx- i 1) r))))))))
  329.  
  330. ;;; Returns a fixnum whose sign is the same as (- |u| |v|).
  331.  
  332. (define (bignum-compare-magnitudes u v)
  333.   (let ((u-size (bignum-length u))
  334.         (v-size (bignum-length v)))
  335.     (cond ((fxn= u-size v-size)
  336.            (fx- u-size v-size))
  337.           (else
  338.            (iterate loop ((i (fx- u-size 1)))
  339.              (receive (diff carry)
  340.                       (%digit-subtract (bignum-digit u i) (bignum-digit v i) 0)
  341.                (cond ((and (fxn= i 0) (fx= diff 0) (fx= carry 0))
  342.                       (loop (fx- i 1)))
  343.                      ((fxn= carry 0) -1)
  344.                      ((fxn= diff 0) 1)
  345.                      (else 0))))))))
  346.  
  347. ;;; Bitfields
  348.  
  349. ;;; only works on positive bignums
  350. (define (bignum-bit-field bn start count)
  351.   (let ((end (fx- (fx+ start count) 1)) )
  352.     (let ((first-hd (fx/ start *bits-per-hyperdigit*))   
  353.           (last-hd  (fx/ end *bits-per-hyperdigit*))
  354.           (bits-in-last (fx+ (fx-rem end *bits-per-hyperdigit*) 1))
  355.           (start-in-hd  (fx-rem start *bits-per-hyperdigit*)) )
  356.       (cond ((or (not (bignum? bn))
  357.                  (not (bignum-positive? bn))
  358.                  (fx< start 0)
  359.                  (fx< count 0)
  360.                  (fx> first-hd (fx+ (bignum-length bn) 1)) 
  361.                  (fx> last-hd  (fx+ (bignum-length bn) 1)) )
  362.              (error "inconsistent arguments in~%  (bignum-bit-field ~s ~s ~s)"
  363.                     bn start count))
  364.             ((fx= first-hd last-hd) 
  365.              (fixnum-bit-field (bignum-digit bn first-hd) start-in-hd count))
  366.             (else
  367.              (let* ((new-last-hd (fx- last-hd first-hd))
  368.                     (new-bn (create-bignum (fx+ new-last-hd 1))))
  369.                (do ((i 0 (fx+ i 1)))
  370.                    ((fx> i new-last-hd)
  371.                     (modify (bignum-digit new-bn new-last-hd)
  372.                             (lambda (hd) (fixnum-bit-field hd 0 bits-in-last)))
  373.                     (set-bignum-sign! new-bn 1)
  374.                     (normalize-integer 
  375.                      (bignum-trim! (bignum-ashr! new-bn start-in-hd)))
  376.                     )
  377.                  (set (bignum-digit new-bn i) (bignum-digit bn (fx+ first-hd i)))
  378.                  )
  379.                ))))))
  380.  
  381.  
  382. ;;; Arithmetic shift entry points:
  383.  
  384. ;;; Left shift bignum.
  385.  
  386. (define (bignum-ashl num amount)
  387.   (receive (q r)
  388.            (%digit-divide 0 amount *bits-per-hyperdigit*)
  389.     (let* ((old-size (bignum-length num))
  390.            (new-size (fx+ (fx+ old-size q) 1)))
  391.       (bignum-trim!
  392.        (bignum-ashl! (make-and-replace-bignum new-size num q 0 old-size)
  393.                      r)))))
  394.  
  395. ;;; Destructive shift.  0 <= amount < *bits-per-hyperdigit*
  396.  
  397. (define (bignum-ashl! num amount)
  398.   (let ((j (fx- (bignum-length num) 1))
  399.         (z (fx- *bits-per-hyperdigit* amount)))
  400.     (iterate loop ((prev (bignum-digit num j))
  401.                    (j j))
  402.       (let ((high (fixnum-ashl prev amount)))
  403.         (cond ((fx<= j 0)
  404.                (set (bignum-digit num j) high)
  405.                num)
  406.               (else
  407.                (let ((foo (bignum-digit num (fx- j 1))))
  408.                  (set (bignum-digit num j) (fx+ high (fixnum-lshr foo z)))
  409.                  (loop foo (fx- j 1)))))))))
  410.  
  411. ;;; Right shift bignum yielding bignum.
  412.  
  413. (define (bignum-ashr src amount)
  414.   (receive (q r)
  415.            (%digit-divide 0 amount *bits-per-hyperdigit*)
  416.     (let* ((old-size (bignum-length src))
  417.            (new-size (fx- old-size q)))
  418.       (bignum-trim!
  419.        (bignum-ashr! (make-and-replace-bignum new-size src 0 q new-size)
  420.                      r)))))
  421.  
  422. ;;; Destructive right shift.
  423.  
  424. (define (bignum-ashr! num amount)
  425.   (let* ((end (fx- (bignum-length num) 1))
  426.          (z (fx- *bits-per-hyperdigit* amount)))
  427.     (iterate loop ((prev (bignum-digit num 0))
  428.                    (j 0))
  429.       (let ((high (fixnum-lshr prev amount)))
  430.         (cond ((fx>= j end)
  431.                (set (bignum-digit num j) high)
  432.                num)
  433.               (else
  434.                (let ((foo (bignum-digit num (fx+ j 1))))
  435.                  (set (bignum-digit num j)
  436.                       (fx+ high (fixnum-ashl foo z)))
  437.                  (loop foo (fx+ j 1)))))))))
  438.